home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_11 / invpow.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  1.4 KB  |  47 lines  |  [MATF/MATL]
  1. function  [lambda,V] = invpow(A,alpha,X,epsilon,max1,show)
  2. % [lambda,V] = invpow(A,alpha,X,epsilon,max1,show)
  3. % The inverse power method is used to find the dominant eigenpair.
  4. % A is an n x n matrix, input.
  5. % X is the starting vector, input.
  6. % alpha is the given shift, input.
  7. % epsilon is the tolerance, input.
  8. % max1 is the maximum number of iterations, input.
  9. % lambda is the dominant eigenvalue, output.
  10. % V is the dominant eigenvector, output.
  11. if nargin==5, show = 0; end
  12. [n n] = size(A);
  13. A = A - alpha*eye(n);               % Form the matrix A - aIpha I
  14. [A,row] = LUfact(A);
  15. lambda = 0;
  16. cnt = 0;
  17. err = 1;
  18. done = 0;
  19. iterating = 1;
  20. state = iterating;
  21. while ((cnt<=max1)&(state==iterating))
  22.   Y = LUsolv(A,X,row);
  23.   [m j] = max(abs(Y));
  24.   c1 = Y(j);                        % Find "largest" element of Y.
  25.   dc = abs(lambda - c1);
  26.   Y = (1/c1)*Y;                     % Perform scalar multiplication
  27.   dv = norm(X-Y);
  28.   err = max(dc,dv);
  29.   X = Y;                            % Update vector  X
  30.   lambda = c1;                      % Update scalar  lambda 
  31.   state = done;
  32.   if (err>epsilon),                 % Check for convergence
  33.     state = iterating;
  34.   end
  35.   cnt = cnt+1;
  36.   lamb = alpha + 1/c1;
  37.   if show==1,
  38.     home; if cnt==1, clc; end; 
  39.     disp(['Inverse Power Method iteration No. ',int2str(cnt)]),...
  40.     disp(''),disp('c1 = '),disp(c1),...
  41.     disp('lambda = '),disp(lamb),...
  42.     disp('Eigenvector = '),disp(X)
  43.   end
  44. end
  45. lambda = alpha + 1/c1;
  46. V = X;
  47.